#include <cmath>
#include "gradient.h"

#include <iostream>
#include <iomanip>
using namespace std;

//функция вычисления градиента
std::valarray<double> gradient (std::valarray<double> w, double (*func)(std::valarray<double>,file_names& Fnames,input& data),
                                input& data, file_names& Fnames) {

    //центральная производная
    unsigned n = w.size();
    std::valarray<double> p(n);
    std::valarray<double> wp(n), wm(n);
    wm = w; wp = w;
    double h0, h, pr0, pr;
    for (unsigned i=0; i<n; i++) {
        h0 = 0.01;
        h = h0;
        wp[i] += h;
        wm[i] -= h;
//            cout << "###" << endl;
//            cout << "f1= " << func(wp,Fnames,data) << endl;
//            cout << "f2= " << func(wm,Fnames,data) << endl;
//            cout << "f1-f2= " << func(wp,Fnames,data)-func(wm,Fnames,data) << endl;
//            cout << "###" << endl;
        pr0 = (func(wp,Fnames,data)-func(wm,Fnames,data))/(2*h0);
        pr = pr0;
        wp[i] = w[i];
        wm[i] = w[i];
        //cout << endl << endl << endl;

        do {
            pr0 = pr;
            h /= 2;
            wp[i] += h;
            wm[i] -= h;
            //cout << "~~~" << endl;
            //cout << "Номер переменной= " << i << endl;
            //cout << "Шаг= " << h << endl;
            //cout << "wp[2]= " << setprecision(15) << wp[2] << endl;
            //cout << "wm[2]= " << setprecision(15) << wm[2] << endl;
            //cout << "f1= " << setprecision(10) << func(wp,Fnames,data) << endl;
            //cout << "f2= " << setprecision(10)<< func(wm,Fnames,data) << endl;
            //cout << "f1-f2= " << setprecision(10) << func(wp,Fnames,data)-func(wm,Fnames,data) << endl;
            //cout << "grad= " << (func(wp,Fnames,data)-func(wm,Fnames,data))/(2*h) << endl;
            //cout << "~~~" << endl;
            pr = (func(wp,Fnames,data)-func(wm,Fnames,data))/(2*h);
            wp[i] = w[i];
            wm[i] = w[i];
        } while (std::abs(pr-pr0) > data.epsilon);
        p[i] = pr;

    }
    return p;
}


/*
    unsigned n = w.size();
    std::valarray<double> p(n);
    std::valarray<double> wh = w;
    float h0, h, pr0, pr;
    for (unsigned i=0; i<n; i++) {
        h0 = 0.01;
        h = h0;
        wh[i] += h;
        pr0 = (func(wh,Fnames,data)-func(w,Fnames,data))/h0;
        pr = pr0;
        wh[i] = w[i];
        do {
            pr0 = pr;
            h /= 2;
            wh[i] += h;
            //for (unsigned i=0; i<n; i++) cout << "wh[" << i << "]= " << wh[i] << endl;
            //cout << "f1= " << func(wh,Fnames,data) << endl;
            //cout << "f2= " << func(w,Fnames,data) << endl;
            pr = (func(wh,Fnames,data)-func(w,Fnames,data))/h;
            wh[i] = w[i];
        } while (std::abs(pr-pr0) > data.epsilon);
        p[i] = pr;

    }
    return p;
*/
